library(dplyr)
library(readxl)
library(regclass)
library(forcats)
library(scales)
library(ggplot2)
library(readxl)
library(writexl)

#The file OE_ReproductionData must first be downloaded - dry season data should then be 
#extracted from sheet 4 and wet season data should be extracted from sheet 5. Data from
#each season should be saved as separate files.
tof.d <- read_excel("REPRODUCTION DATA LOCATION - DRY SEASON")
tof.w <- read_excel("REPRODUCTION DATA LOCATION - WET SEASON")

####Variable preparation####

#dry season
tof.d$north.south <- as.factor(tof.d$north.south)
tof.d$use.category2 <- as.factor(tof.d$use.category2)
tof.d$mango.yn <- as.factor(tof.d$mango.yn)
tof.d$bg.yn <- as.factor(tof.d$bg.yn)
tof.d$ban.yn <- as.factor(tof.d$ban.yn)
tof.d$pap.yn <- as.factor(tof.d$pap.yn)
tof.d$avo.yn <- as.factor(tof.d$avo.yn)
tof.d$ca.yn <- as.factor(tof.d$ca.yn)
tof.d$guava.yn <- as.factor(tof.d$guava.yn)
tof.d$gmel.yn <- as.factor(tof.d$gmel.yn)


tof.d$edu.lvl <- factor(tof.d$edu.lvl, levels = c("None", "Primary", "Secondary", "Diploma"))
tof.d$edu.lvl <- fct_collapse(tof.d$edu.lvl, Secondary = c("Secondary", "Diploma"))

tof.d$hh.size_s <- rescale(tof.d$hh.size)
tof.d$crop.count_s <- rescale(tof.d$crop.count)
tof.d$no.acres.cult_s<- rescale(tof.d$no.acres.cult)
tof.d$TLU_s<- rescale(tof.d$TLU)
tof.d$hh.mpi_s <- rescale(tof.d$hh.mpi)
tof.d$FC.1km_s <- rescale(tof.d$FC.1km)
tof.d$kcal.EUI.F_s <- rescale(tof.d$kcal.EUI.D)

#wet season

tof.w$north.south <- as.factor(tof.w$north.south)
tof.w$use.category2 <- as.factor(tof.w$use.category2)
tof.w$mango.yn <- as.factor(tof.w$mango.yn)
tof.w$bg.yn <- as.factor(tof.w$bg.yn)
tof.w$ban.yn <- as.factor(tof.w$ban.yn)
tof.w$pap.yn <- as.factor(tof.w$pap.yn)
tof.w$avo.yn <- as.factor(tof.w$avo.yn)
tof.w$ca.yn <- as.factor(tof.w$ca.yn)
tof.w$guava.yn <- as.factor(tof.w$guava.yn)
tof.w$gmel.yn <- as.factor(tof.w$gmel.yn)

tof.w$edu.lvl <- factor(tof.w$edu.lvl, levels = c("None", "Primary", "Secondary", "Diploma"))
tof.w$edu.lvl <- fct_collapse(tof.w$edu.lvl, Secondary = c("Secondary", "Diploma"))

tof.w$hh.size_s <- rescale(tof.w$hh.size)
tof.w$crop.count_s <- rescale(tof.w$crop.count)
tof.w$no.acres.cult_s<- rescale(tof.w$no.acres.cult)
tof.w$TLU_s<- rescale(tof.w$TLU)
tof.w$hh.mpi_s <- rescale(tof.w$hh.mpi)
tof.w$FC.1km_s <- rescale(tof.w$FC.1km)
tof.w$kcal.EUI.W_s <- rescale(tof.w$kcal.EUI.W)

#creating the use category factor variable, using "No/Other trees" as the reference level

tof.d <- within(tof.d, use.category2 <- relevel(use.category2, ref = "NOT"))
tof.w <- within(tof.w, use.category2 <- relevel(use.category2, ref = "NOT"))

####NUTRITION MODELS####

##MEASURING ASSOCIATIONS BETWEEN ON-FARM TREE USE AND MICRONUTRIENT ADEQUACY
#code provided for first dry and then wet season datasets

tof.z.d <- lm(Zinc.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(tof.z.d)

tof.z.w <- lm(Zinc.MSM.NAR ~ use.category2 + FC.1km_s  + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(tof.z.w)

tof.va.d <- lm(VitA.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(tof.va.d)

tof.va.w <- lm(VitA.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(tof.va.w)

tof.i.d <- lm(Iron.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(tof.i.d)

tof.i.w <- lm(Iron.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(tof.i.w)

tof.f.d <- lm(Folate.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(tof.f.d)

tof.f.w <- lm(Folate.MSM.NAR ~ use.category2 + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(tof.f.w)


####FOOD SOURCING MODELS####

##MEASURING ASSOCIATIONS BETWEEN ON-FARM TREE USE AND FOOD SOURCING PATTERNS
#code provided for first dry and then wet season datasets

wild.tof.d <- lm(wild.pct ~ use.category2 + kcal.EUI.D_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(wild.tof.d)

wild.tof.w <- lm(wild.pct ~ use.category2 + kcal.EUI.W_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(wild.tof.w)

cult.tof.d <- lm(cult.all.pct ~ use.category2 + kcal.EUI.D_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(cult.tof.d)

cult.tof.w <- lm(cult.all.pct ~ use.category2 + kcal.EUI.W_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(cult.tof.w)

pur.tof.d <- lm(purchased.pct ~ use.category2 + kcal.EUI.D_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(pur.tof.d)

pur.tof.w <- lm(purchased.pct ~ use.category2 + kcal.EUI.W_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(pur.tof.w)


####TABLE S7####
#Linear models testing the effect of cultivating different species of on-farm trees on 
#women’s micronutrient adequacy (vitamin A, zinc, iron, folate) in both dry and wet seasons.

#NOTE: each model can be run for each micronutrient outcome (all are currently measuring Zinc)

#mango
test.mango.d <- lm(Zinc.MSM.NAR ~ mango.yn + tree.count_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.mango.d)

test.mango.w <- lm(Zinc.MSM.NAR ~ mango.yn + tree.count_s + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.mango.w)

#eucalyptus ("bluegum")
test.bg.d <- lm(Zinc.MSM.NAR ~ bg.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.bg.d)

test.bg.w <- lm(Zinc.MSM.NAR ~ bg.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.bg.w)

#banana
test.ban.d <- lm(Zinc.MSM.NAR ~ ban.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.ban.d)

test.ban.w <- lm(Zinc.MSM.NAR ~ ban.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.ban.w)

#papaya
test.pap.d <- lm(Zinc.MSM.NAR ~ pap.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.pap.d)

test.pap.w <- lm(Zinc.MSM.NAR ~ pap.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.pap.w)

#gmelina
test.gmel.d <- lm(Zinc.MSM.NAR ~ gmel.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.gmel.d)

test.gmel.w <- lm(Zinc.MSM.NAR ~ gmel.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.gmel.w)

#avocado
test.avo.d <- lm(Zinc.MSM.NAR ~ avo.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.avo.d)

test.avo.w <- lm(Zinc.MSM.NAR ~ avo.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.avo.w)

#custard apple
test.ca.d <- lm(Zinc.MSM.NAR ~ ca.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.ca.d)

test.ca.w <- lm(Zinc.MSM.NAR ~ ca.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.ca.w)

#guava
test.guava.d <- lm(Zinc.MSM.NAR ~ guava.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.d)
summary(test.guava.d)

test.guava.w <- lm(Zinc.MSM.NAR ~ guava.yn + FC.1km_s + north.south + hh.size_s + hh.mpi_s + edu.lvl + no.acres.cult_s + crop.count_s + TLU_s, data = tof.w)
summary(test.guava.w)

